<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Figures 6.21-6.23: Basis pursuit using Gabor functions</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/cvxbook/Ch06_approx_fitting/html/basispursuit.html">
<link rel="stylesheet" href="../../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Figures 6.21-6.23: Basis pursuit using Gabor functions</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Section 6.5.4</span>
<span class="comment">% Boyd &amp; Vandenberghe "Convex Optimization"</span>
<span class="comment">% Original by Lieven Vandenberghe</span>
<span class="comment">% Adapted for CVX by Argyris Zymnis - 11/27/2005</span>
<span class="comment">%</span>
<span class="comment">% Here we find a sparse basis for a signal y out of</span>
<span class="comment">% a set of Gabor functions. We do this by solving</span>
<span class="comment">%       minimize  ||A*x-y||_2 + ||x||_1</span>
<span class="comment">%</span>
<span class="comment">% where the columns of A are sampled Gabor functions.</span>
<span class="comment">% We then fix the sparsity pattern obtained and solve</span>
<span class="comment">%       minimize  ||A*x-y||_2</span>
<span class="comment">%</span>
<span class="comment">% NOTE: The file takes a while to run</span>

clear

<span class="comment">% Problem parameters</span>
sigma = 0.05;  <span class="comment">% Size of Gaussian function</span>
Tinv  = 500;   <span class="comment">% Inverse of sample time</span>
Thr   = 0.001; <span class="comment">% Basis signal threshold</span>
kmax  = 30;    <span class="comment">% Number of signals are 2*kmax+1</span>
w0    = 5;     <span class="comment">% Base frequency (w0 * kmax should be 150 for good results)</span>

<span class="comment">% Build sine/cosine basis</span>
fprintf(1,<span class="string">'Building dictionary matrix...'</span>);
<span class="comment">% Gaussian kernels</span>
TK = (Tinv+1)*(2*kmax+1);
t  = (0:Tinv)'/Tinv;
A  = exp(-t.^2/(sigma^2));
ns = nnz(A&gt;=Thr)-1;
A  = A([ns+1:-1:1,2:ns+1],:);
ii = (0:2*ns)';
jj = ones(2*ns+1,1)*(1:Tinv+1);
oT = ones(1,Tinv+1);
A  = sparse(ii(:,oT)+jj,jj,A(:,oT));
A  = A(ns+1:ns+Tinv+1,:);
<span class="comment">% Sine/Cosine basis</span>
k  = [ 0, reshape( [ 1 : kmax ; 1 : kmax ], 1, 2 * kmax ) ];
p  = zeros(1,2*kmax+1); p(3:2:end) = -pi/2;
SC = cos(w0*t*k+ones(Tinv+1,1)*p);
<span class="comment">% Multiply</span>
ii = 1:numel(SC);
jj = rem(ii-1,Tinv+1)+1;
A  = sparse(ii,jj,SC(:)) * A;
A  = reshape(A,Tinv+1,(Tinv+1)*(2*kmax+1));
fprintf(1,<span class="string">'done.\n'</span>);

<span class="comment">% Construct example signal</span>
a = 0.5*sin(t*11)+1;
theta = sin(5*t)*30;
b = a.*sin(theta);

<span class="comment">% Solve the Basis Pursuit problem</span>
disp(<span class="string">'Solving Basis Pursuit problem...'</span>);
tic
cvx_begin
    variable <span class="string">x(30561)</span>
    minimize(sum_square(A*x-b)+norm(x,1))
cvx_end
disp(<span class="string">'done'</span>);
toc

<span class="comment">% Reoptimize problem over nonzero coefficients</span>
p = find(abs(x) &gt; 1e-5);
A2 = A(:,p);
x2 = A2 \ b;

<span class="comment">% Constants</span>
M = 61; <span class="comment">% Number of different Basis signals</span>
sk = 250; <span class="comment">% Index of s = 0.5</span>

<span class="comment">% Plot example basis functions;</span>
<span class="comment">%if (0) % to do this, re-run basispursuit.m to create A</span>
figure(1); clf;
subplot(3,1,1); plot(t,A(:,M*sk+1)); axis([0 1 -1 1]);
title(<span class="string">'Basis function 1'</span>);
subplot(3,1,2); plot(t,A(:,M*sk+31)); axis([0 1 -1 1]);
title(<span class="string">'Basis function 2'</span>);
subplot(3,1,3); plot(t,A(:,M*sk+61)); axis([0 1 -1 1]);
title(<span class="string">'Basis function 3'</span>);
<span class="comment">%print -deps bp-dict_helv.eps</span>

<span class="comment">% Plot reconstructed signal</span>
figure(2); clf;
subplot(2,1,1);
plot(t,A2*x2,<span class="string">'--'</span>,t,b,<span class="string">'-'</span>); axis([0 1 -1.5 1.5]);
xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'y_{hat} and y'</span>);
title(<span class="string">'Original and Reconstructed signals'</span>)
subplot(2,1,2);
plot(t,A2*x2-b); axis([0 1 -0.06 0.06]);
title(<span class="string">'Reconstruction error'</span>)
xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'y - y_{hat}'</span>);
<span class="comment">%print -deps bp-approx_helv.eps</span>

<span class="comment">% Plot frequency plot</span>
figure(3); clf;

subplot(2,1,1);
plot(t,b); xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'y'</span>); axis([0 1 -1.5 1.5]);
title(<span class="string">'Original Signal'</span>)
subplot(2,1,2);
plot(t,150*abs(cos(w0*t)),<span class="string">'--'</span>);
hold <span class="string">on</span>;
<span class="keyword">for</span> k = 1:length(t);
  <span class="keyword">if</span>(abs(x((k-1)*M+1)) &gt; 1e-5), plot(t(k),0,<span class="string">'o'</span>); <span class="keyword">end</span>;
  <span class="keyword">for</span> j = 2:2:kmax*2
    <span class="keyword">if</span>((abs(x((k-1)*M+j)) &gt; 1e-5) | (abs(x((k-1)*M+j+1)) &gt; 1e-5)),
      plot(t(k),w0*j/2,<span class="string">'o'</span>);
    <span class="keyword">end</span>;
  <span class="keyword">end</span>;
<span class="keyword">end</span>;
xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'w'</span>);
title(<span class="string">'Instantaneous frequency'</span>)
hold <span class="string">off</span>;
</pre>
<a id="output"></a>
<pre class="codeoutput">
Building dictionary matrix...done.
Solving Basis Pursuit problem...
 
Calling sedumi: 61625 variables, 502 equality constraints
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 502, order n = 61125, dim = 61626, blocks = 2
nnz(A) = 7484105 + 0, nnz(ADA) = 252004, nnz(L) = 126253
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            7.63E+04 0.000
  1 :  -1.00E-01 3.44E+02 0.000 0.0045 0.9990 0.9990   3.00  1  1  5.1E-01
  2 :   5.06E+00 1.08E+02 0.000 0.3136 0.9000 0.9000   1.00  1  1  4.8E-01
  3 :   1.03E+01 4.82E+01 0.000 0.4468 0.9000 0.9000   1.00  1  1  4.1E-01
  4 :   1.19E+01 2.73E+01 0.000 0.5657 0.9000 0.9000   1.00  1  1  3.5E-01
  5 :   1.25E+01 1.52E+01 0.000 0.5567 0.9000 0.9000   1.00  1  1  2.7E-01
  6 :   1.27E+01 9.41E+00 0.000 0.6201 0.9000 0.9034   1.00  1  1  2.1E-01
  7 :   1.27E+01 5.70E+00 0.000 0.6057 0.9000 0.9068   1.00  1  1  1.5E-01
  8 :   1.28E+01 2.83E+00 0.000 0.4968 0.9272 0.9000   1.00  1  1  9.1E-02
  9 :   1.28E+01 1.33E+00 0.000 0.4707 0.9496 0.9000   1.00  1  1  4.7E-02
 10 :   1.28E+01 9.07E-01 0.000 0.6808 0.9000 0.9081   1.00  1  1  3.3E-02
 11 :   1.28E+01 5.50E-01 0.000 0.6068 0.9325 0.9000   1.00  1  1  2.1E-02
 12 :   1.28E+01 3.38E-01 0.000 0.6146 0.9170 0.9000   1.00  1  1  1.3E-02
 13 :   1.28E+01 1.99E-01 0.000 0.5880 0.9353 0.9000   1.00  1  1  7.6E-03
 14 :   1.28E+01 1.19E-01 0.000 0.6001 0.9052 0.9000   1.00  1  1  4.6E-03
 15 :   1.28E+01 7.74E-02 0.000 0.6485 0.9037 0.9000   1.00  1  1  3.0E-03
 16 :   1.28E+01 3.41E-02 0.000 0.4402 0.9000 0.9042   1.00  1  1  1.3E-03
 17 :   1.28E+01 1.26E-02 0.000 0.3690 0.9184 0.9000   1.00  2  2  4.9E-04
 18 :   1.28E+01 1.83E-03 0.000 0.1454 0.9190 0.9000   1.00  2  2  7.1E-05
 19 :   1.28E+01 4.35E-04 0.000 0.2377 0.9000 0.9070   1.00  2  2  1.7E-05
 20 :   1.28E+01 8.31E-05 0.000 0.1912 0.9011 0.9000   1.00  2  2  3.2E-06
 21 :   1.28E+01 2.17E-05 0.000 0.2607 0.9000 0.9018   1.00  2  2  8.5E-07
 22 :   1.28E+01 7.19E-06 0.000 0.3320 0.9000 0.9031   1.00  2  2  2.8E-07
 23 :   1.28E+01 2.79E-06 0.000 0.3884 0.9006 0.9000   1.00  3  3  1.1E-07
 24 :   1.28E+01 1.07E-06 0.000 0.3825 0.9079 0.9000   1.00  3  3  4.2E-08
 25 :   1.28E+01 3.84E-07 0.000 0.3589 0.9112 0.9000   1.00  3  3  1.5E-08
 26 :   1.28E+01 1.38E-07 0.000 0.3596 0.9076 0.9000   1.00  4  4  5.4E-09

iter seconds digits       c*x               b*y
 26     94.0   8.0  1.2845155824e+01  1.2845155686e+01
|Ax-b| =   1.6e-11, [Ay-c]_+ =   1.1E-11, |x|=  2.6e+00, |y|=  9.9e-01

Detailed timing (sec)
   Pre          IPM          Post
2.820E+00    9.395E+01    3.000E-01    
Max-norms: ||b||=1.497811e+00, ||c|| = 1,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 3306.09.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +12.8452
done
Elapsed time is 86.389939 seconds.
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="basispursuit__01.png" alt=""> <img src="basispursuit__02.png" alt=""> <img src="basispursuit__03.png" alt=""> 
</div>
</div>
</body>
</html>